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This document presents the results of the research work performed by Drs. G. S. Liaw and 
J. D. Mo of Alabama A&M university for NASA/Marshall Space Flight Center, under 
Grant NAG8-064. Two major tasks were accomplished during the course of this research 
project. First, a new computer code was developed for the low thrust viscous nozzle 
flowfield predictions by using a new LU scheme. This task was proposed at the beginning 
of this contract Second, the existing FDNS code was implemented to include the radiation 
effect. This task was added into this project later. Results for a Carbon Dioxide and the 
SSME nozzles are documented in this final report. 
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I. Introduction 


This research project has three objectives. The first objective is to compute the Carbon 
Dioxide nozzle flow fields by using an existing computer code. This work has been 

completed and documented previouslyfl]. The second objective is to develop a new 
computer code for solving the compressible full Navier-Stokes equations for low thrust 
type nozzle flow calculations. The low density effect is embedded through the no-slip 
conditions on the wall boundaries. The third objective is to analyze the space shuttle Main 
Engine (SSME) flow field in the combustor and the divergent nozzle with radiative effect. 
The radiation model has been derived and embedded into the existing FDNS computer code 
which is currently operational in ED32 NASA/MSFC. However, our main contribution is 
to develop a compressible Navier-Stokes code which is capable of handling the rarefication 
effect 

Basically, there are two kind of numerical schemes to solve the time-dependent fluid 
dynamics problems. They are explicit or implicit. The early developments, limited by the 
computer capability, emphasized on the explicit methods. Several explicit numerical 
schemes have been developed successfully, such as the Euler explicit method, the Lax 
method, Leap frog method, MacCormack method, etc.. The obvious way to accelerate 
convergence to a steady state or to save the computer time for unsteady problems is to 
increase the size of the time step. However, it was found that the time step size for the 
explicit schemes is seriously limited by the Courant-Friedrichs-Lewy (CFL) condition, 
which requires that the region of dependence of the difference scheme must be a subset of 
the region of dependence of the differential equation. 

Generally speaking, implicit schemes are preferred when the time step limit imposed by an 
explicit bound is much less than that imposed by the accuracy bound. The computation is 
unconditionally stable for implicit schemes, and the time step is determined by the desired 
level of accuracy. However, the implicit methods require to solve a large number of 
coupled equations at each time step. Hence, the reduction in the number of time steps may 
be outweighed by the increase in the number of arithmetic operations required for each time 
step. As the computer technology advances, the restrictions on the storage and computing 
time have been relaxed. The implicit methods become more popular in the CFD 
community. 

In this work, a Navier-Stokes code has been developed for the low thrust viscous nozzle 
flow fields prediction. An implicit finite volume, in an arbitrary curvilinear coordinate 
system, lower-upper(LU) scheme[2] is used to solve the governing Navier-Stokes 
equations and species transportation equations. This scheme was originally developed by 
Jameson, and extended to an axisymmetric coordinate system by this group. Sample 
calculations of Carbon Dioxide nozzle flow are presented in this report to verify the validity 
and efficiency of this code. The computed results are in reasonable agreement with the 
experimental data. The bench marie date were chosen from Chou and Carter[3]. 
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II. The LU-Code Development 


2.1. Formulation 

Basic governing equations: T he governing Navier-Stokes equations are derived from the 
basic physical laws of conservation of mass, momentum and energy. These equations are 
cast into the conservative forms in the cylindrical coordinate system 
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This set of compressible flow equations govern the steady or unsteady; laminar or 
turbulent, chemically reacting or nonreacting flow problems. S is the summation of all the 

source terms, while 5=0 for the rectangular coordinate and 5=1 for the axisymmetric 
coordinate system. 
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Grid generation and coordinate transformation: The numerical solution of the system of the 
governing partial differential equations can be greatly simplified by a well constructed grid 
network. It is also true that a grid which is not well suited to the problem can produce 
unsatisfactory results. Improper choice of grid point locations may lead to an apparent 
instability. 

Early numerical solution by finite difference methods was restricted to problems where 
suitable coordinate systems must be selected in order to solve the governing equations in 
this compatible coordinate system. As more complex flowfields problems are under 
consideration, general mappings have been employed to transform the physical plane into a 
computational domain. Numerous advantages accrue when appropriate transfer procedures 
are followed, for example, the body surface can be selected as a boundary in the 
computational plane permitting easy application of surface boundary conditions. Also, in 
general, transformation is used which lead to a uniformly spaced grid in the computational 
plane while points in physical space may be unequally spaced. 

For general purposes, consider the following general transformation, which can be 
orthogonal or nonorthogonal, 

£ =£ (x,y) (22) 

ri=ri(x, y) (23) 

the computational grid points are located along the transformed coordinates (E,, rj), which 
can be designed to have both equal or unequal spacing in the computational plane 

depending on the problems. 

Generally, the grid generation can be accomplished by three different approaches. 

* the complex variable method, 

* the algebraic method and 

* the differential equations method. 

The complex variable technique has the advantage that the transformations used are analytic 
or partially analytic as opposed to those methods that are entirely numerical. 
Unfortunately, complex variable methods are restricted in two dimensional problem. For 
this reason, the technique has limited applicability and will not be used in the present work. 
Algebraic and differential equation techniques can be used on both two and three 
dimensional problems with complex geometry, they are adopted in this code development. 

To be consistent with the general grid system, a corresponding coordinate transformation 
on the governing partial differential equation is required. The requirements of this 
transformation for grid generation are: 

* The mapping is one - one correspondent, i.e. no singularity exists in the 
mapping function. 

* The grid lines is smooth to assure that the derivatives of the mapping function are 
continuous. 

* Grid points are closely clustered in the region where the large numerical errors 
are expected. 

* Avoid excessive grid skewness. 

For the transformation (2.2) and (2.3), the governing equations are converted from the 
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physical domain (x, y) to the computational domain ( q). By using the chain rule of 
partial differentiation, the partial derivative becomes 
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( 2 . 4 ) 
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The matrices ( 2; Xj £y, q Xf rjy ) appearing in these equations can be determined in the 
following manner: 
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In a similar manner, we can write 
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where the Jacobian 
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The terms x^, x^ , yj-, and y^ are obtained directly by the central difference. Then, the 

quantities % x , q X) r|y, which appear in the governing equation, are evaluated from 
equation (2.10). 
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After the coordinate transformation, the governing equation (2.1) becomes: 
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The strong conservation form of the above equations are 
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It can be written in a compact form. 
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where 


( 2 . 12 ) 


f=| x F+| y G, f v =^ x F v+ ^ y G v 

g —T\ x F +7) y G, fy =7J x Fy + t) yGy 

and S the source term is unchanged by the transformation. 

2.2. Numerical Procedures 

Linearization of the governing equations bv the delta form formulation: For simplicity, the 
derivation of the equations below is carried out in the Cartesian coordinate system, then the 
results are transformed onto a general coordinate system. 
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Using the implicit scheme with central differencing, Eq. (2.12) can be formulated as 


Q n+1 = Q n -/3 <5 1 jr?c (f( Q n+1 )-f v ( Q n+1 )+<?», (g( Q n+1 ) — g v < Q n+1 )-S n+1 ] 

-(i-fhStfc (ft Q n )-f v (Q ")+<?,, (g(Q n )-g v <Q n )-S n | (2.13) 


where 8 £ and 8^ are difference operators of 8/8!;, 8/8r|, and P is a positive number 

between 0 and 1, which is used as an adjust parameter, with P=0 for an explicit scheme 

and P=1 for an implicit scheme. P=1/2 is designated for the Crank-Nicolson scheme 
which has a second order accuracy in the time discretization. 

The Jacobian matrices are 
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and the increment of the conservative variable Q is 

5 Q =Q n+1 -Q n (2.15) 

The scheme is linearized by setting 

f(Q n+1 )=f(Q n ) Q + 0(|£ Q f ) 

g(Q n+1 )=g(Q")+|^5Q+0(|«Q| 2 ) 

a y 

S(Q n+1 )=S(Q n )+^-|5Q+0(|5Q| 2 ) (2-16) 

a y 


where terms of the second and higher order have been omitted. This yields, 
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In compact form, the final 5-form of equation (2.17) becomes 

jl+ /? At (<9c AA +<?,, BB-h)}+ At R=0 (2.18) 


where 
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( 2 . 19 ) 
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scheme. We begin with the following procedures, by setting. 
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(2.23) 
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(2.24) 

with the definitions of matrices 


AA+= 1/2(AA+v a I) 

(2.25) 

AA- = 1/2 (AA -v A I) 

(2.26) 

then 


9^ AA = 9^ + AA' + 9^' AA+ 

(2.27) 

In the similar manner, 


9 n AA = 9 r j + BB- + 9rj' BB+ 

(2.28) 

where. 


BB+ = 1/2 ( BB + v B I ) 

(2.29) 

BB' = 1/2 (BB - v B I ) 

(2.30) 

After substitution, equation (3.7) becomes 


jl + 0 At (<9f AA + +<?^AA " + + + 9* BB” -H )j SQ =-AtR 

(2.31) 

where 9p’ and 9^’ are defined as backward-difference operators and 9 f + and 

defined as forward difference operators. 


The values v A and v B are chosen to construct AA + , AA', 

BB + , BB' so 


eigenvalues of "+" matrices are nonnegative and those of-"-" matrices are nonpositive. 
The development of these matrices is extremely important for the success of the LU-type 
scheme. 
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Jameson and Yoon [4] chose 


v A > max ( |X A |) (2.32) 

v B > max ( I X B | ) (2.33) 


where X A and Xg are the eigenvalues of matrix AA and BB. 
After some manipulations, equation (2.31 ) becomes: 
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(2.34) 


From (2.25), (2.26), (2.29) and (2.30), 


AA+-AA' = va I 
BB+-BB* = v0| 


Equation (2.34 ) can be factored in the following forms 

jl + )3 At[(AA + — AA “)+(BB + -BB - )j+/? At |^'AA + + ^"BB + -AA + -BB + -H jj 

*|l + )3 At|(AA + -AA ")+(BB + -BB')]+)3 Al(^ + AA'+^ + BB - +AA'+BB')J<5Q 
= -|l+ )3At(u A +u B )]AtR (2.35) 


This factorization is arbitrary as long as the final numerical scheme is stable. On the other 
hand, a further improvement on the scheme can be done from this point. 

Numerical discretization and computational procedures: Based on the discretized governing 
equation (2.35), a finite-volume method is applied to discrete the spatial variable and 
separated time and spatial discretization is involved to assure a steady state solution 
independent of the time step. The finite volume formulation provides a convenient treatment 
of complex geometry and avoids the problems of metric singularity which are usually 
associated with the finite difference methods. 

The control volume is shown in the Fig. 1. Point P is chosen as the control point located at 
the centroid of the control volume. 
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Fig. 1. Control volume and control point 

By integrating the equation (2.35) over the control volume, the left hand side of the 
equation take the value at the control point as its average value. The first two terms in the 
right hand side constitute a divergent vectors and they are envoluated by a surface 
integration after applying the Gauss theorem. The residual source term takes the value at the 
control point. 

In order to match the approximated factorization of the lower-upper triangular matrix form, 
the derivatives in the left hand side is discretized by the first order backward and forward 
difference schemes, such that 




(2.36) 
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where the subscripts (i,j) represent the location of the grid point (^,rj) and superscript 
stands for the time level, in which superscript n means at the time t and n+1 means at time 

t+ At. is the projected spatial distance between the control point (i- 1 ,j) and (i,j) in £ 
coordinate direction and Arj is the projected spatial distance between the control point (i,j-l) 
and (i,j) in rj coordinate direction. A£ + and Ar| + are similarly defined. 

After the discretization of the partial differential equations, they become a system of 
algebraic equations. The continuous governing partial differential equations have been 
replaced by a set of discretized finite difference equations which are valid only at the 
discretized grid points. At each grid point, a finite difference equation is created and the 
connection among all the grid points in the solution domain is hold by the forward and 
backward finite difference operators. With proper boundary conditions at the boundaries 
of the solution domain and initial conditions at the beginning of the physical process for the 
unsteady problems, the set of the algebraic equations can be solved by proper numerical 
procedures. 

A mathematical problem is described by the following algebraic equations to illustrate the 
essence of the lower-upper factorization procedures of the matrix conversion. 

A X = b (2.38) 

If the coefficient matrix A can be decomposed into a lower-upper matrix multiplication 
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form, then the equation (2.38) can be rewritten as 

LU X = b (2.39) 

Then we can solve this equation in two steps. At first, let Y=U X, equation (3.31) 
becomes 


L Y = b (2.40) 

Because the matrix L is a lower triangular matrix. Its inverse L'* can be easily determined 
by the forward Gaussian elimination method, or 

Y = L-lb (2.41) 

Then, looking for the solution for x such that 

UX=Y (2.42) 

Since U is an upper triangular matrix, the matrix inverse can be found by the one step 
backward Gaussian elimination method. 


X = U* 1 Y 


(2.43) 


The lower-upper factorization procedure in the present research follows the same 
procedures as described above. The dimension of the final coefficient matrix is 


[ C ] = [ (I x J) x (IxJ) ] 

where I and J are the grid number in the % and r| directions respectively. 

Boundary conditions: Two types of geometrical boundaries are considered. One is solid 
walls, the other is non-solid walls which includes inlets, the symmetric axis, and planes of 
symmetry. 

Solid walls: The no-slip conditions are applied on any solid walls, . In addition, zero 
pressure gradient are employed. 


Inlet: For generality, any part of the boundaries is referred as inlet where the flow is 
inward to the solution domain. At the inlet, the flow variables are assumed known. That 
is, the velocities, static pressure and temperature are specified. 

Exit: In the same way, any part of the boundary is classified as exit where the flow is 
outward to the solution domain. For generality, the exit is divided into two categories, 
which are supersonic or subsonic exits. For supersonic exits, no boundary conditions are 
needed because its hyperbolic nature. For subsonic exits, back pressure must be specified, 
and other flow variables apply 



(2.44) 


where <|) represents any flow variable except the static pressure. 

It should be noted that even for supersonic exit, a boundary layer exists near the exit wall. 
So a known backpressure is applied in the subsonic boundary layer region. 
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Initial condition: Our interest is to achieve steady slate solutions even though the 
procedures described above are valid for general unsteady flow problems. Therefore, only 
an appropriate initial solution is needed to start the time-iterative scheme, and finally the 
steady state solution is reached through the time marching process. For nozzle flows, an 
one dimensional inviscid supersonic solution is used as the initial condition. 

2.3. Computer Program 

The computer program has been structured in a general manner. At the present state, the 
program can treat compressible flows, turbulent as well as lamilar, chemically reactive as 
well as non-reactive, internal as well as external, two-dimensional or axisymmetric 
arbitrary geometry problems. Variable transport properties are also included. 

The program structure is written in a modular form so that subroutines for specific physical 
effects can be flagged only when those effects are needed. After its execution, the main 
processor writes a considerable amount of information to external data files. This 
information is used by the postprocessor to produce a variety of computer graphic displays. 
These data will be made use of to restart the calculations in the future work. Most physical 
and numerical control parameters can be changed which gives the users considerable 
flexibility for their specific purposes, the basic flow diagram of the computer program is 
shown in Fig. 2. 

2.4. Carbon Dioxide Nozzle Row Calculation 

For nozzle flowfield calculations, five reservoir pressures (29.4, 14.7, 7.35, 3.7,1.85)psi 
were run and shown in Fig.3. The reservoir temperature is fixed at 600°F. Due to the 
cryogenic flow conditions in downstream, all the variable transport properties including the 
specific heat ratio are curve-fitted and extrapolated to extremely low temperature, shown in 
Fig.4. As for the boundary layer development at different chamber pressures, the velocity 
profiles and the boundary thicknesses at the nozzle exit are shown in Fig.5(a) and 5(b), 
respectively. The boundary layer thicknesses is defined at the boundary edge where the 
velocity magnitude is 99% of the core flow. It is found that the lower the chamber 
pressure, the thicker the boundary layer at the nozzle exit. This flow characteristics is due 
to the Reynolds number resulted from the low chamber pressure.For the chamber pressure 
1.85psi, the boundary layer thickness is about one third of nozzle exit radius. As a result, 
the conventional MOC/BL concept is no longer applicable. The wall pressures of both 
computed and test data are shown in Fig.6. Fig. 6(a) and 6(b) show the wall pressure 
distributions in linear and logarithmic scales. They are in good agreement. Test data 
asymptotically deteriorate as the exit pressures and temperatures approach the liquid-vapor 
saturation line in the CC> 2 -phase diagram. Condensation phenomena in the supersonic gas 
stream remains an unexplored research topics. For chamber pressure less than 1.85, both 
rarefication effect and cryogenic effect create physical as well as numerical instability. Our 
computer experiments to impose an explicit slip wall boundary conditions have failed. The 
implicit boundary conditions become the necessary step for next trials. 


III. FDNS Code Implementation 


Currently, the FDNS Code has been comprehensively used for the prediction of the SSME 
nozzle characteristics in NASA/MSFC[5]. The basic equations employed in the code are 
the axisymmetric, multi-component conservation equations. A generalized form of these 
equations written in curvilinear coordinates is given by 
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where J, Uj and Gj j represent the Jacobian of the coordinate transformation, contravariant 
velocities and diffusion metrics respectively. 

An unpaid scheme was employed to approximate the convective terms of the momentum, 
energy and continuity equations. The scheme is based on second and fourth order central 
differencing with artificial dissipation. First order upwinding is used for the species and 
turbulence equations, since the parameters involved are positive quantities. Different 
eigenvalues are used for weighing the dissipation terms depending on the conserved 
quantity being evaluated, in order to give correct diffusion fluxes near wall boundaries. 

For the SSME combuster, the temperature is in the value of 6000°R. It is natural to ask if 
the wall heat flux includes the radiation contribution. We have added this capability into the 
FDNS code and delivered it to MSFC ED32, and results were published in ref.[6]. 

The radiation heat transfer has been hypothesized by a six-flux radiation model, which is 
based on the Schuster-Hanaker approximation in astrophysical research. For an 
absorbing-emitting gray medium in local thermodynamic equilibrium, the radiation 
transport equations describing the variations of the fluxes along six direction can be 
reduced to the following three second-order ordinary differential equations: 
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Where the composite-fluxes R x , Ry, and Rq are defined as : 

R *=j(V +I r) 

R r4(V +I r) 


In these equations, I x +, I r +, Iq+, I x -, I r and Iq are the radiation energy fluxes along 

the positive and negative coordinate directions. The scattering coefficient, s, is defined as 
the scattering radiation energy per unit length. 


The FDNS has been coded and implemented. The calculations have been conducted for the 
SSME combustor and nozzle. Fig. 7 shows the configuration and the computational grid. 
The grid was generated with more grid point near the wall to predict the large gradient 
there. Fig. 8 shows the isotherm and velocity vector plot of the entire flowfield and Fig.9 
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shows the averaged radiation flux contours and velocity vector plot. Fig. 10 shows the 
enlarged averaged radiation flux contours in the throat area. Comparisons have been made 
for the heat flux over the wall and temperature along the centerline with and without 
radiation contribution. The results are shown in Fig. 1 1 and Fig. 12, respectively. 


IV. Personnel 

In this project, Dr. G. S. Liaw serves as the Principal Investigator to coordinate all the 
research activities. Dr. J. D. Mo, a Research Assistant Professor, has been working in full 
time to develop the LU-code since September 15, 1989 until the project was Finished. 


V. Conclusions 

A new computer code for analyzing the axisymmetric nozzle flow with variable gamma has 
been developed. The validity of this code is demonstrated by the comparisons of present 
calculations with experimental data. The code has been used to simulate the flowfield in a 
Carbon Dioxide nozzle having an area ratio of 40. The radiation modeling within the SSME 
combustor/nozzle has shown that the radiation heat transfer has relatively significant 
contribution along the nozzle throat, but little effect downstream. The total effect on the 
wall heat flux is within 5%, which is consistent with the qualitative analysis from MSFC in 
house results. The code has shown to be efficient and accurate for the flow conditions 
considered in the present study. 
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Fig. 2. Row chart of LU procedure 
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chamber pressure: i) 294, ii) 14 7. iii) 7.35, 
iv) 3 7 and v) 1.85 psi 


ii U1 IV 

normalized velocity profiles at exit 


chamber pressure, psi 


Fig. 5. Velocity profiles and boundary layer thickness 
at nozzle exit 
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Fig. 6. Comparison of computed wall pressure 
and test data 


Fig. 7. SSME combustor/nozzle configuration and 
computational grid 



Fig. 8. Isothermal and velocity vector plot of flowfield 
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Fig. 9. Averaged radiation flux contours and 
velocity vector plot 



Fig. 10. Enlarged averaged radiation flux contours and 
velocity vector plot in the throat area 
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Fig. 11. The relative difference of heat flux over the wall 
with and without radiation contribution 



Fig. 12. The relative difference of the temperature 
along the centerline 
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